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ABSTRACT 


This  paper  presents  the  methodology  and  results  of  a computational  experiment 
which  compares  the  performance  of  four  computer  codes  which  determine  the  best 
discrete  approximation  to  a continuous  nonlinear  function.  The  experiment 
utilizes  320  test  problems  created  by  a test  problem  generator.  Several 
performance  measures  describe  solution  quality  as  well  as  computational 
effort. 


Key  Words:  Algorithm  testing;  approximation;  computational  experiment;  least 

absolute  deviation;  performance  measures;  polynomial  approximation. 
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I.  INTRODUCTION 

This  paper  is  part  of  a series  of  papers  dealing  with  computational  experi- 
mentation in  mathematical  programming,  which  emphasizes  sound,  statistically- 
based  methodology  for  comparing  codes  in  this  area.  As  the  specific  object  of 
our  experimentation  (and,  therefore,  a vehicle  for  developing  a sound  method- 
ology), we  chose  to  study  codes  that  solve  problems  of  least  absolute  value 
curve  fitting.  Estimation  in  the  norm  provides  a useful  adjunct  to  least 
squares  estimation  [12,20],  especially  when  errors  derive  from  long-tailed 
distributions  [15,16]. 

In  earlier  work  [14],  we  studied  four  codes,  representative  of  a range  of  so- 
lution techniques,  using  a statistically  based  experiment  design  and  pseudo- 
randomly  generated  test  problems  (representative  of  general  data-fitting 
problems).  It  was  discovered  that  such  codes  could  be  clearly  ranked  with  re- 
spect to  computational  effort  (CPU  time),  but  that  very  little  difference 
among  the  codes  could  be  observed  with  respect  to  certain  "quality-of- 
results"  performance  indicators:  accuracy,  consistency,  and  correctness. 

Therefore,  another  experiment  was  designed  using  polynomial  approximation  test 
problems  and  the  same  four  codes  used  in  the  earlier  work.  Polynomial  ap- 
proximation problems  are  frequently  encountered  in  practice  [7,19,25],  and  are 
known  to  admit  some  of  the  interesting  properties  of  ill-conditioning  [11] 
which  tend  to  differentiate  among  codes  with  respect  to  the  quality  of  re- 


sults 
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Moreover,  prior  to  executing  the  first  comparison  experiment,  we  tested  the 
codes  on  a battery  of  hand-picked  test  problems  obtained  from  the  authors,  the 
literature  and  code  users  [22].  Each  of  the  codes  yielded  incorrect  informa- 
tion on  at  least  one  of  the  hand-picked  problems,  with  polynomial  approxima- 
tion problems  being  especially  troublesome. 

This  paper  will  present  an  experiment  for  testing  codes  on  a class  of  poly- 
nomial approximation  problems,  and  discuss  the  results  of  this  test  effort. 
Directions  for  future  research  in  both  code  development  and  code-comparison 
methodology  will  be  presented  in  the  final  section  of  this  report. 
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II.  BACKGROUND 


Given  n sets  of  observations  (yi,x-j_)  i=l,  • . . , n,  the  discrete  L]  polynomi- 
al approximation  problem  is  that  of  determining  values  0O,  • . . , 0m  which 
minimize 


n m j 

0 (8q»  • • • » Bm)  3 I I yi  — L x 

i-1  j-0  1 


(1) 


In  other  words,  for  n observations  on  a dependent  variable  y and  an  explana- 
tory variable  x,  we  wish  to  determine  a polynomial  of  degree  m which  best  fits 
these  data,  in  the  sense  of  minimizing  the  sum  of  the  absolute  values  of  the 
residuals 


j 

ei  “ 7i  ” I 8j  x 
j=0 


if  it 

A vector  6*  = (0  , . . . , 0 ) which  yields  the  minimum  value  0*  of 
o m 

0 (0O,  . . . , 0m)  is  termed  a solution  vector,  with  optimum  objective  func- 


tion value  0*.  Unlike  the  case  for  polynomial  approximation  over  a contin- 


uous interval  [25],  page  38,  the  solution  vector  for  a discrete  polynomial 


approximation  problem  need  not  be  unique. 


Discrete  polynomial  approximation  problems  can  be  formulated  and  solved  as 


linear  programming  problems  by  the  stratagem  of  introducing  variables  a 

i 
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minimize 


n 


I 

i=l 


ai 


m 3 

ubject  to  x > Yi  , 


=0  J i 


i = 1 , . . . , n , 


(2) 


m 


ai  l 0j 
j=0  i 


> - 


Yi 


i — 1 , . . . , n • 


It  can  be  demonstrated  [2,4,13]  that  a solution  vector  3*  can  always  be  found 
that  interpolates  at  least  m+1  of  the  data  points;  that  is,  at  least  m+1  re- 
siduals e-j_  are  equal  to  zero  at  the  solution.  Degeneracy  is  said  to  occur 
when  more  than  m+1  residuals  have  value  zero  at  a solution  B*. 


Certain  algorithms  for  solving  general  discrete  approximation  problems  make 
use  of  the  above  "interpolation"  result.  Namely,  these  algorithms  restrict 
attention  to  basic  solutions  (where  precisely  m+1  residuals  equal  zero),  and 
move  in  a systematic  way  from  one  basic  solution  to  an  improved  basic  solu- 
tion. Any  such  collection  of  m+1  zero-valued  residuals  corresponds,  then,  to 
the  specification  of  m+1  active  constraints:  values  selected  from  i 3 1 , . . 


, n,  where  both  inequalities  of  (2)  hold  with  equality. 
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Four  general-purpose  discrete  approximation  FORTRAN  codes  were  selected  for 
study  and  evaluation  in  this  paper.  Three  of  these  codes  are  based  on  linear 
programming  reformulations  of  the  approximation  problem: 

(I)  A double-precision  dual  revised  simplex  code  by  Abdelmalek  [1], 

(II)  A primal  revised  simplex  code  with  a partial  sort  procedure  by  Arm- 
strong [3], 

(III)  A primal  simplex  code  with  full  tableau  and  limited  use  of  double- 
precision by  Barrodale  and  Roberts  [5]. 

While  codes  II  and  III  are  essentially  primal  simplex  codes,  they  do  make  use 
of  certain  dual  properties  in  passing  over  several  simplex  vertices  in  a sin- 
gle iteration. 

In  addition,  a restricted  gradient  method  for  minimizing  piecewise  differenti- 
able functions  is  embodied  in  a fourth  code: 

(IV)  A descent  code  by  Bartels  [6],  incorporating  a form  of  QR  decom- 
position. 

A May  1977  version  of  each  code  was  obtained  from  the  authors.  Tolerance  set- 
tings were  fixed  at  those  levels  suggested  by  the  authors  as  appropriate  for 
our  computer  (UNIVAC  1108  with  36-bit  words).  In  code  IV,  a random  perturba- 
tion is  internally  applied  to  the  problem  data  when  near-degeneracy  is  suspected. 
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Unfortunately,  when  using  this  code  as  a subroutine  to  solve  a series  of  identical 
problems,  the  actual  perturbation  applied  varies  from  problem  to  problem.  In  order 
to  ensure  that  the  same  (perturbed)  problem  is  solved  in  successive  executions 
using  the  same  input  problem,  the  pseudo-random  number  generator  provided  in 
Code  IV  was  always  initialized  prior  to  execution.  Other  than  these  changes 
(tolerance  settings,  and  pseudo-random  number  initialization),  the  codes  were 
implanted  exactly  as  provided  by  the  authors.  Certain  general  characteristics 
of  the  various  codes  are  indicated  in  Table  I. 


Table  I 


General  Code  Characteristics 


Code 

Code 

Code 

Code 

I 

II 

III 

IV 

Number  of 
Statements 

270 

186 

223 

1609 

Required* 

Storage 

2np+2p2+7n 
+6p+l 15 1 

np+p^+6n 

+4p+824 

np+5p+5n 

+723 

np+1 . 5p2+5n 
+6. 5p+4496 

Restrictions 

None 

Full  Rank 

None 

tn>  1 

Additional 

Features 

None 

None 

None 

Allows  constraints 
variable  starting 
points,  selective 
printing 

Defined  in 

terms  of  n = number 

of  observations 

, and  p =■ 

m+1  * number  of 

estimated  parameters. 
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In  our  previous  study  of  these  codes  [14],  two  types  of  test  problems  were 
considered:  (1)  a collection  of  27  hand-picked  problems  [22],  devised  with  a 

specified  structure  in  mind  or  arising  from  an  actual  application;  and  (2) 
pseudo-randomly  generated  problems  for  general  discrete  approximation  [17]. 
While  for  certain  classes  of  pseudo-randomly  generated  problems,  all  four 
codes  "correctly"  solved  each  generated  problem,  this  was  not  the  case  for 
many  of  the  27  hand-picked  problems.  Since  a number  of  these  "difficult" 
problems  were  in  fact  forms  of  discrete  polynomial  approximation  in  the 
norm,  it  appeared  useful  to  study  systematically  the  performance  of  the  four 
codes  on  a wider  class  of  polynomial  approximation  problems.  To  this  end, 
test  problem  generators  have  been  developed  [9]  to  create  a range  of  polynomi- 
al approximation  problems  representative  of  situations  encountered  in  prac- 
tice. 

© 

One  test  problem  generator  (P0LY1)  models  the  situation  in  which  a given  non- 
linear function  f(x)  is  to  be  approximated  in  the  norm  by  a polynomial  of 
specified  degree,  over  a discrete  set.  Accordingly,  the  observed  values  y^  in 
(2)  are  modeled  as  arising  from  the  given  function,  via  yj_  = f(x^).  Control- 
lable features  of  this  generator  include:  the  choice  of  function  f(x),  the 

interval  I of  approximation,  the  distribution  and  number  of  observation  points 
Xj_el,  and  the  degree  of  the  approximating  polynomial.  For  this  particular 
generator,  the  actual  solution  vector  8*  cannot  be  specified  in  advance.  How- 
ever, in  the  second  generator  (P0LY2)  for  another  class  of  polynomial  approxi- 
mation problems,  the  solution  vector  can  be  so  specified;  the  current  experi- 
ments have  concentrated  on  the  use  of  P0LY1,  and  thus  attention  will  subse- 
quently be  confined  to  results  obtained  using  this  first  generator. 
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With  P0LY1  one  can  study  discrete  polynomial  approximation  of  a variety  of 
functional  forms,  defined  on  certain  intervals.  The  eight  functions  f(x)  em- 
ployed in  the  present  code  evaluations  are  listed  in  Table  II,  together  with 
their  associated  domains.  These  functions,  and  their  intervals  of  definition, 
were  chosen  in  order  to  represent  a useful  variety  of  shapes  likely  to  be  en- 
countered and  (through  the  choice  of  interval)  a range  of  "ill-conditioning” 
for  the  problem  matrix  X = (x^).  Figure  1 graphs  the  light  functions  studied 
in  this  paper. 


TABLE  II 

Functions  and  Their  Intervals  of  Approximation 
Function  f(x)  Interval  I 


1 

e“x  sin  x 

. [0,4] 

2 

ex  sin  x 

[0,4] 

3 

ex  sin  x 

[1,7] 

4 

e2x/ 2x 

[-05,1 

5 

75x/[l+(7.5x)2] 

[0,2] 

6 

10  x e“*^x 

[0,4] 

7 

l/[l+(x-2.5)4] 

[0,5] 

8 

xl/3 

[-1,1] 

9 


» » S 3 « 

X 


2.  Y=eXsin(x) 


5.  Y=75x[l+(7-5x)2]  6<  y=lOxe  * 5x 
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III.  METHODOLOGY 

The  approach  taken  in  comparing  the  four  codes  consisted  of  designing  a 
computational  experiment  to  discern  the  effects  of  certain  controlled  vari- 
ables by  examining  a number  of  performance  indicators.  This  section  of  the 
paper  will  describe  the  details  of  this  experiment. 

The  first  step  in  designing  the  experiment  involved  the  identification,  speci 
fication,  and  quantification  of  the  controlled  variables  which  characterize 
the  problems  generated  by  P0LY1.  These  controlled  variables  fall  into  three 
main  categories:  those  describing  the  function  to  be  approximated,  those  de- 

scribing the  spacing  of  the  observations  within  the  interval  over  which  che 
function  is  approximated,  and  those  describing  the  size  of  the  problem.  For 
the  experiment  described  here,  the  functions  used  were  the  eight  functions 
described  in  Table  II  together  with  their  associated  intervals  of  approxima- 
tion. Observations  were  chosen  to  be  evenly  spaced  within  the  interval. 

The  size  of  each  approximation  problem  was  specified  by  the  number  of  observa 
tions  n and  the  degree  m of  the  approximating  polynomial.  For  the  experiment 
the  number  of  observations  assumed  the  values  50,  100,  200,  400,  and  800  and 
the  approximating  polynomial  had  degree  varying  from  3 to  10  inclusive.  Thus 
each  of  the  four  codes  was  exercised  on  a total  of  320  problems  (8  func- 
tions, 8 degrees,  5 sets  of  equally  spaced  observations). 
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The  second  part  of  the  experiment  design  involved  identifying  and  quantifying 
various  indicators  of  code  performance.  Performance  indicators  can  be  grouped 
naturally  into  two  general  categories:  those  related  to  computational  effort 

and  those  related  to  the  quality  of  the  solution  produced.  The  indicators  in 
both  categories  must  be  judged  relative  to  the  particular  computer  system  and 
computer  environment  in  which  the  experiment  is  conducted.  All  of  the  results 
reported  in  this  paper  were  obtained  on  the  UNIVAC  1108  at  the  National  Bureau 
of  Standards  under  the  Exec  8 operating  system,  Version  33. R3,  with  36  bits 
for  each  single  precision  computer  word. 

Perhaps  the  most  frequently  used  performance  indicator  is  Central  Processing 
Unit  (CPU)  time.  In  this  experiment  CPU  time  was  chosen  as  the  appropriate 
indicator  of  computational  effort  inasmuch  as  the  four  codes  were  sufficiently 
different  for  iteration  counts  to  be  meaningless.  Quantification  of  this  in- 
dicator was  obtained  by  accessing  an  internal  clock  from  a system  subroutine 
which  references  an  appropriate  elapsed  CPU  time  entry  (in  milliseconds)  in 
the  run  accounting  table.  Previous  results  in  timing  experiments  on  our  com- 
puter [13]  indicated  several  sources  of  variability  in  the  process  which  con- 
found code  comparison  results.  First,  the  inherent  variability  associated 
with  clock  and  accounting  table  update  was  found  to  be  on  the  order  of  2 or  3 
milliseconds.  A second  possible  contributor  to  variability  occurs  because  the 
execution  time  depends  on  the  locations  in  core  memory  of  instructions  and  the 
data  upon  which  those  instructions  operate.  To  reduce  this  effect,  the  dimen- 
sions of  the  arrays  used  by  all  four  codes  are  fixed  and  do  not  vary  with 
problem  size.  A third  source  of  timing  variability  is  the  computing  environ- 
ment: the  number  and  types  of  other  jobs  being  run  during  the  same  time 
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period.  Previous  results  [14]  indicated  that  for  any  environment  other  than  a 
totally  dedicated  computer,  one  run  of  a code  on  a problem  would  not  suffice 
to  determine  true  CPU  time  on  that  problem.  Thus,  all  of  the  timing  results 
reported  here  were  obtained  under  totally  dedicated  conditions. 

The  second  category  of  performance  indicators  studied  relates  to  the  quality 
of  the  solution  obtained.  This  category  (which  is  not  always  considered) 
proved  to  be  critically  important  in  the  present  experiment,  whereas  in  our 
previous  experiment  [14],  it  did  not  play  a decisive  role  in  distinguishing 
among  the  codes.  The  evaluation  of  solution  quality  requires  either  knowledge 
of  the  "true"  solution  or  the  ability  to  check  a given  solution  for  optimali- 
ty. As  mentioned  in  Section  2,  the  polynomial  approximation  problems  produced 
by  P0JLY1  cannot  be  generated  with  solutions  known  a priori.  Thus,  it  was 
necessary  to  devise  a procedure  for  checking  solution  quality.  Accordingly, 
we  developed  a double  precision  computer  subroutine  which,  given  the  problem 
data  and  specified  set  of  active  constraints,  checks  the  Kuhn-Tucker  condi- 
tions for  optimality  [4,17].  If  satisfied,  it  then  performs  a double- 
precision reinversion  of  the  optimal  basis  matrix  to  obtain  accurate  optimal 
values  for  the  objective  function  and  the  solution  vector.  These  solution 
values  are  accepted  as  the  "true"  values  for  purposes  of  comparison  with  the 
results  produced  by  the  codes. 

Four  different  performance  indicators  were  chosen  to  represent  solution  quali- 
ty: correctness,  accuracy,  accuracy  after  reinversion,  and  consistency.  \1 1 

four  indicators  involve  comparisons  between  "true"  and  code-derived  active 
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constraints,  solution  vectors  or  objective  function  values.  Several  quantifi- 
cations were  used  in  this  study.  For  comparing  objective  function  values,  the 
experiment  calculated  the  number  of  digits  of  agreement,  the  absolute  differ- 
ence, and  the  relative  difference.  For  comparing  solution  vectors,  the  ex- 
periment examined  average  number  of  digits  of  agreement,  absolute  Euclidean 

* 

distance,  and  normalized  Euclidean  distance  (where  the  two  vectors  Bco<je  and 
* 

Btrue  were  viewed  as  points  in  (m-i-1) -dimensional  space).  Precise  definitions 
of  these  indicators  and  their  quantification  are  deferred  to  the  next  section. 

The  experiment  design  reported  here  thus  consists  of  the  specification  of 
problems  to  be  solved  by  each  of  the  four  codes  and  the  identification  of 
the  indicators  upon  which  code  performance  will  be  evaluated  and  compared. 

The  next  section  presents  the  results  of  executing  this  experiment  design. 
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IV.  RESULTS 

Each  of  the  four  codes  was  tested  using  all  320  test  problems.  Elapsed  CPU 
time  and  four  measures  for  quality  of  results  (correctness,  accuracy,  accuracy 
after  reinversion,  and  consistency)  were  obtained. 


CORRECTNESS 

Correctness  (did  the  code  actually  arrive  at  the  true  solution?)  is  a measure 
of  performance  of  interest  to  all  users.  We  defined  correctness  as  a binary 
(YES/NO)  measure.  If  the  active  constraints  defining  the  final  solution  pro- 
duced by  an  Lj_  code  satisfied  the  Kuhn-Tucker  conditions  to  given  tolerances, 
it  is  called  "correct".  In  our  analyses,  we  noticed  that  some  codes  produced 
an  objective  function  value* which  agreed  with  the  true  solution  to  eight  stg- 
nificant  digits  but  whose  final  simplex  vertex  was  not  an  optimal  vertex. 

Since  comparisons  using  more  than  eight  significant  digits  are  inappropriate 
using  single-precision  arithmetic  on  the  UNIVAC  L108,  whenever  the  objective 
function  produced  by  a single  precision  code  agreed  with  that  of  the  true 
solution  to  eight  digits,  we  also  considered  the  solution  "correct”.  For  ten 
of  the  320  problems,  no  code  produced  the  correct  solution  judged  according  to 
Kuhn-Tucker  conditions.  These  problems  were  eliminated  from  our  subsequent 
analysis  since  the  "true"  solution  could  not  be  ascertained.  Analyses  of  the 
results  using  the  correctness  measure  are  presented  in  Figures  2-4.  There 
is  no  table  for  Code  I which  obtained  the  correct  solution  for  every  problem 
analyzed  (310  of  the  320).  Code  II  produced  incorrect  results — regardless  of 
function  type — for  almost  all  problems  of  degree  greater  than  5 or  for 
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Figure  2. 
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Figure  3. 

CODE  III  -CORRECT  SOLUTIONS 
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Figure  '4. 
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problems  with  greater  than  100  observations.  Codes  III  and  IV  were  less  fre- 
quently correct  than  Code  I but  more  often  correct  than  Code  II.  Although  one 
can  conclude  that  both  of  these  codes  are  likely  to  be  judged  "not  correct" 
for  problems  having  both  high  degree  and  many  observations,  consistent  pat- 
terns did  not  emerge  to  delineate  the  function-type  or  problem  sizes  which  are 
most  likely  to  cause  difficulty  to  either  of  these  codes. 


ACCURACY 


Having  established  that  all  codes  do  not  yield  the  correct  solution  to  all 
problems,  the  next  natural  question  is  how  close  to  correctness  are  the  an- 
swers provided  by  the  codes.  Accordingly,  the  performance  indicator,  accur- 
acy, compares  the  true  solution  with  that  obtained  by  the  code.  Figures 
5-7  present  for  each  code  the  fraction  of  solutions  produced  with  errors  In* 
the  objective  function  values  no  greater  than  a given  magnitude,  where  error 
is  quantified  by  three  measures:  the  absolute  difference,  the  relative  dif- 

ference, and  the  number  of  digits  of  agreement  between  the  true  objective 
function  value  and  that  produced  by  the  code.  The  test  results  for  all  310 
test  problems  produces  a definitive  ranking  of  the  codes,  with  Code  I being 
the  superior  code,  followed  by  Code  III,  then  Code  IV  and  finally  Code  II. 

Note  that  this  ranking  is  invariant  with  respect  to  the  measurement  quantifi- 
cation. 

Three  measures  were  developed  for  quantifying  the  accuracy  of  the  solution 
vector.  The  true  solution  vector  and  that  produced  by  each  code  can  be 
considered  as  points  in  (m+l)-space.  The  Euclidean  distance  and  the 


Figure  5. 
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normalized  Euclidean  distance  between  the  points  are  then  two  measures  of  the 
accuracy  of  the  solution  vectors  produced  by  the  codes.  A third  measure  is 
obtained  by  calculating  the  number  of  digits  of  agreement  (0  to  8)  between 
corresponding  components  of  the  solution  vectors  and  averaging  these  over  Che 
m+1  components  to  obtain  the  average  number  of  accurate  digits  in  the  solution 
vector.  These  three  measures  of  solution  vector  accuracy  are  presented  in 
Figures  8 - 10.  Note  that  the  ranking  of  the  four  codes  is  essentially  the 
same  as  that  obtained  for  objective  function  accuracy  except  that  Codes  III 
and  IV  are  more  evenly  matched. 

ACCURACY  AFTER  REINVERSION 


In  performing  the  analyses  of  these  test  results,  we  noticed  that  frequently  a 
code  would  arrive  at  a "correct"  vertex  but  the-  objective  function  value  and/ 
or  the  solution  vector  would  be  "inaccurate".  We,  therefore,  performed  a 
double-precision  inversion  of  the  final  basis  matrix  and  compared  the  objec- 
tive function  value  and  solution  vector  obtained  after  inversion  to  the  true 
solution  values.  We  call  this  performance  measure  "accuracy  after  reinver- 
sion" . 

We  note  that  for  Code  IV  the  results  here  are  based  on  only  279  (rather  than 
310)  problems  because  31  of  the  solutions  produced  by  Code  IV  had  fewer  than 
m+1  active  constraints,  thereby  rendering  reinversion  impossible.  Figure  11 
presents  the  results  of  examining  the  relative  differences  in  objective  func- 
tion values.  Similar  results  were  obtained  for  the  other  two  quantifications, 
i.  e.  absolute  difference  and  number  of  digits  of  agreement.  Figure  12 
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Figure  9. 
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illustrates  the  accuracy  in  the  solution  vector  after  reinversion  as  quanti- 
fied by  the  normalized  Euclidean  distance  between  the  two  points.  Again,  sim- 
ilar results  were  obtained  for  the  absolute  Euclidean  distance  and  the  average 
number  of  digits  of  agreement.  Note  that  the  accuracy  of  codes  II,  III,  and 
IV  have  improved  substantially  through  reinversion.  Thus,  if  one  were  to  use 
any  of  the  single  precision  codes  tested  in  the  study,  it  is  clear  that  an  ac- 
curate (double  precision)  inversion  of  the  final  basis  is  desirable. 

CONSISTENCY 


A fourth  indicator  of  solution  quality,  consistency , provides  some  information 
about  the  effects  of  round-off  error.  It  refers  to  the  internal  consistency 
between  the  solution  vector  and  the  objective  function  value  produced  by  the 
code.  The  solution  vector  produced  by  the  code  is  used  with  the  given  problem 
data  to  calculate,  in  double-precision,  the  residuals  and  an  objective  func- 
tion value.  Consistency  is  then  measured  by  comparing  this  objective  function 
value  with  that  given  by  the  code.  As  can  be  seen  from  Figure  13,  approxi- 
mately 40  percent  of  the  solutions  produced  by  Codes  II,  III,  and  IV  yielded 
inconsistencies  in  the  solution  information  such  that  the  relative  differences 
in  objective  function  value  exceeded  10“^.  The  ranking  of  the  codes  for  this 
performance  measure  is  not  unlike  that  produced  for  the  other  performance  mea- 
sures. Note  that  even  the  "correct"  and  "accurate"  results  of  Code  I are  not 
consistent,  an  indication  that  round-off  error  is  affecting  the  results. 


Figure  13. 

RELATIUE  CONSISTENCY  OBJECTIUE  FUNCTION  UAIUE 


29 


apruruSexi  usaiq 

uet[i  on  Joaj;a  suoxunxog  jo  ijuaoxaa 


30 


In  summary,  one  observes  that  Code  I — a double-precision  code — was  superior  to 
the  other  three  single-precision  codes  with  respect  to  all  measures  of  all 
four  quality  of  results  criteria.  It  should  also  be  noted  that  Code  III, 
which  is  a single-precision  code  using  one  double-precision  variable  to  ac- 
cumulate inner  products,  generally  outperforms  the  other  two  single-precision 
codes.  This  inexpensive  and  core-conserving  device  should  be  considered  in 
all  single-precision  codes  where  numerical  difficulties  may  be  anticipated. 

Finally,  although  two  of  the  three  single-precision  codes  had  a switch  to  In- 
dicate that  they  terminated  due  to  numerical  difficulties  and  although  there 
were  many  instances  when  significant  round-off  errors  existed,  only  rarely  did 
a code  terminate  with  this  warning. 

TIMING 

As  stated  in  Section  IIT,  CPU  time  for  runs  made  on  a dedicated  computer  was 
used  as  an  indicator  of  computational  effort.  In  a previous  comparison  ex- 
periment [14]  with  the  four  codes,  solution  quality  was  "perfect"  and  the 
only  mechanism  for  distinguishing  among  the  codes  was  in  terms  of  CPU  time. 

The  results  of  that  experiment  indicated  that  Code  II  always  required  the 
least  amount  of  time,  followed  by  Codes  III,  I,  and  IV,  generally  in  that 
order.  For  the  experiment  described  herein,  the  solution  quality  became  of 
significant  importance  and,  as  a result,  CPU  time  became  less  important  since 
it  can  be  argued  that  solution  time  is  irrelevant  when  solution  quality  is  van- 
acceptable.  It  can  also  be  argued,  however,  that  CPU  time  is  always  important 
since  (1)  it  is  perhaps  the  single  most  important  factor  In  calculating 
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computer  run  costs  under  many  computer  accounting  systems,  (2)  the  codes  rare- 
ly produced  warnings  that  solution  quality  was  perhaps  suspect,  and  (3)  one 
does  not  know  in  advance  of  code  execution  just  what  to  expect  in  terms  of  the 
solution  time  and  quality.  Thus,  CPU  times  were  analyzed  for  three  of  the 
four  codes  studied  in  this  experiment.  Code  II  was  omitted  from  the  analysis 
of  CPU  time  because  of  the  very  large  number  of  problems  for  which  solution 
quality  was  deemed  unacceptable. 

CPU  times  for  Codes  I,  III,  and  IV  were  analyzed  for  all  320  problems  (in- 
cludes the  10  problems  for  which  no  code  found  a correct  solution)  with  one 
general  exception.  There  were  a total  of  18  problems  for  which  Code  IV  ap- 
parently entered  some  sort  of  infinite  loop,  stopping  the  solution  process  on- 
ly by  exceeding  a maximum  of  1000  iterations  without  determining  a result 
which  the  code  considered  a solution.  These  eight  problems  required  a total 
of  slightly  over  13  minutes  of  computer  time  and  were  not  included  in  the  re- 
sults shown  for  Code  IV.  Because  Code  IV  generally  requires  more  CPU  time 
than  the  other  codes,  the  omission  of  these  exceptionally  large  times  does  not 
significantly  alter  the  ranking  of  the  codes  with  respect  to  CPU  time. 

Figures  14,  15,  and  16  illustrate  various  results  with  respect  to  CPU  time. 
Figure  14  illustrates  the  average  CPU  time  in  milliseconds  as  it  relates  to 
the  number  of  observations.  The  results  here  indicate  that  with  the  possible 
exception  of  a very  slight  reversal  of  the  order  between  Codes  I and  II  at  50 
observations.  Code  I requires  the  least  CPU  time  followed  by  Codes  III  and  IV. 


Figure  14. 

AVERAGE  CPU  TIME  PER  PROBLEM  VS.  NO.  OF  OBSERVATIONS 
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AVERAGE  CPU  TIME  PER  PROBLEM  VS.  FUNCTION 
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In  Figure  15,  the  results  are  not  quite  so  definitive.  This  figure  illus- 
trates average  CPU  time  as  it  relates  to  the  degree  of  the  approximating  poly- 
nomial. Here,  Code  IV  performs  well  for  polynomials  of  low  degree,  but  uses 
noticeably  more  time  for  degrees  9 and  10.  In  contrast,  Codes  I and  III  per- 
form fairly  well  and  consistently  for  all  degrees  studied.  Except  for  degree 
3,  the  average  CPU  time  for  Code  I was  always  less  than  that  for  Code  III. 

Figure  16  illustrates  the  relationship  between  average  CPU  time  and  the  func- 
tion to  be  approximated.  This  figure  illustrates  quite  clearly  the  need  for 
comparing  code  performance  over  a wide  class  of  problems.  Consider,  for  exam- 
ple, the  difference  of  opinion  that  would  result  if  one  were  to  examine  code 
performance  only  on  function  3 or  only  on  function  5 versus  the  examination  of 
all  8 functions.  Indeed,  to  the  extent  that  longer  CPU  time  is  indicative  of 
greater  problem  "difficulty”,  the  results  here  demonstrate  conclusively  that 
problems  which  are  the  most  difficult  for  one  code  are  not  necessarily  the 
most  difficult  (and  at  times  are  even  the  least  difficult)  for  some  other 
code.  Thus,  valid  statements  regarding  comparative  code  performance  cannot  be 
made  solely  on  the  basis  of  computational  experience  with  a small  number  of 
problems. 
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V.  OTHER  RESULTS 

A preliminary  investigation  was  also  made  of  the  effect  of  "ill-conditioning" 
on  the  nature  of  test  problems  and  the  behavior  of  L[  codes  on  those  problems. 
For  the  purpose  of  this  investigation,  the  1-condition  number  (or  H-condition 
number)  of  the  optimal  basis  matrix  B was  chosen  as  a measure  of  problem 
difficulty  [10,18].  An  estimate  k(B)  of  this  condition  number  can  be  effi- 
ciently computed  [10],  with  the  interpretation  that  calculation  of  6*,  using 
Bg*  = YB»  involves  the  loss  of  approximately  logigk(B)  significant  figures  of 
accuracy. 

It  was  found  that  the  condition  number  estimate  k(B)  did  not  vary  appreciably 

with  n (the  number  of  observations),  but  it  was  strongly  affected  by  m (the 

degree  of  the  approximating  polynomial).  In  fact,  for  each  unit  increment  in 

m,  k(B)  grows  by  approximately  one  order  of  magnitude.  Furthermore,  for  fixed 

n and  m,  the  condition  estimate  varies  by  function  - with  function  3 giving 

the  largest  and  function  8 the  smallest  such  estimates.  This  occurence  Ls  not 

unexpected,  in  that  function  3 is  approximated  over  a fairly  broad  interval 

[0,7],  while  function  8 is  approximated  over  a fairly  narrow  interval  [-1,1]. 

More  generally,  k(B)  was  observed  to  increase  with  the  increasing  magnitude  of 

the  interval  endpoints.  Such  results  are  consistent  with  the  notion  that  the 

m 

successive  powers  1,  x^,  . . . , x^  comprising  matrix  B become  more  widely 
divergent  in  magnitude  as  m increases  or  Ix-jJ  increases. 

In  order  to  assess  how  the  accuracy  of  the  codes  changed  In  response  to 

problem  characteristics,  the  average  number  of  agreeing  digits  (AAD)  between 

* * 
the  true  solution  vector  Btrue  anc*  provided  by  the  code  3C0(je  was 
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computed.  Agreement  was  invariably  perfect  for  Code  I,  and  so  subsequent 
analysis  focused  on  Codes  II  - IV.  For  these  latter  codes,  AAD  did  not  vary 
in  a discernable  way  with  n,  but  generally  exhibited  a decrease  with  increas- 
ing values  of  m.  While  this  behavior  of  AAD  parallels  that  found  for  k(B), 
there  was  no  consistent  relation  between  this  measure  of  accuracy  and  the  con- 
dition estimate  for  a given  degree  m.  In  other  words,  very  little  additional 
information  on  code  accuracy  is  conveyed  by  knowing  the  condition  estimate,  in 
addition  to  the  value  of  m. 

REGULARLY  SPACED  POINTS  OF  EXACT  FIT 


The  test  problems  used  within  this  study  were  generated  by  varying  the  form  of 
the  function  to  be  approximated,  the  degree  of  the  approximating  polynomial, 
and  the  number  of  observation  points.  The  observation  points  were  positioned 
evenly  spaced  over  the  interval  of  interest. 

When  the  interpolation  points  defining  the  best  L^  solution  to  each  problem 
were  tabulated — keeping  the  degree  of  the  polynomial  and  the  number  of  obser- 
vaton  points  fixed  but  allowing  the  function-type  to  vary — a pattern  emerged 
which  indicated  that  these  interpolation  points  would  be  regularly  spaced 
within  the  interval  regardless  of  the  functional  form  chosen. 

An  understanding  of  approximation  theory  presents  a partial  explanation  for 
this  phenomenon.  A theorem  in  this  field  [21,  pp.  71-73]  states  that  the  best 
Li  approximating  polynomial  of  degree  m-1  (or  less)  to  x^  over  [-1,1]  has  its 
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interpolation  points  at  the  values  x^  = cos  k tt  (k  = 1,  . . . , m).  These 

m+1 

values  x^  are  the  roots  of  a Chebyshev  polynomial  of  the  second  kind. * This 
theorem  only  holds  for  continuous  approximation  and  describes  Che  fitting  of  a 
polynomial  of  degree  m-1  to  a polynomial  of  next  highest  degree. 

In  our  experiment,  we  were  using  only  a discrete  number  of  points  (evenly 
spaced)  and  non-polynomial  functions.  However,  our  preliminary  analyses  have 
indicated  that  although  the  roots  of  the  Chebyshev  polynomial  are  not  always 
the  interpolation  points  defining  the  best  solution,  they  are  usually  very 
"close".  This  result  suggests  that  one  might  consider  using  such  values  x^  as 
starting  points  for  discrete  polynomial  approximation  algorithms.  Further 
research  in  the  area  is  presented  in  [8]. 
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A mth  order  Chebyshev  polynomial  of  the  second  kind  has  the  following  form: 


m/2  (1-m)! 

Um(x)  = l (-1)1 (2x)ra-2l 

1=0  1 ! (m-21) ! 


where  Um(l)  is  defined  to  equal  m+1. 
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VI.  SUMMARY  AND  RECOMMENDATIONS 

In  summarizing  the  results  of  the  work  described  here,  several  conclusions  and 
recommendations  can  be  made.  First,  the  results  clearly  indicate  that  solu- 
tion quality  is  substantially  improved  by  a final  (double-precision)  reinver- 
sion of  the  optimal  basis  matrix.  We  recommend  that  this  procedure  be  imple- 
mented to  supplement  any  of  the  three  single  precision  Lj.  codes  when  solving 
polynomial  approximation  problems  of  the  type  described  here. 

In  an  earlier  experiment  [14]  using  Lj_  problems  with  different  structure,  all 
four  codes  consistently  produced  correct  solutions  so  that  CPU  time  was  the 
only  distinguishing  performance  measure.  Code  II  used  significantly  less  CPU 
time  than  any  of  the  other  three  codes.  In  particular,  Code  I,  the  double- 
precision code,  was  relatively  slow  in  obtaining  solutions.  The  results  from 
the  polynomial  approximation  experiment  contrast  sharply.  Solution  quality  is 
not  consistent  among  codes.  In  fact,  based  on  the  measures  of  solution  quali- 
ty presented  in  this  report,  it  can  be  concluded  that  Code  II  is  not  well- 
suited  for  solving  the  type  of  polynomial  approximation  problems  utilized  in 
this  experiment.  On  the  other  hand,  Code  I not  only  produced  more  correct 
solutions  than  any  of  the  other  three  codes,  but  it  utilized  relatively  little 
CPU  time  in  doing  so.  It  is  conjectured  that  the  use  of  double-precision 
arithmetic  in  Code  I was  costly  in  terms  of  increased  CPU  time  in  the  first 
experiment,  but  paid  off  in  the  second  experiment  where  the  generated  problems 
were  numerically  more  difficult  to  solve. 
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Codes  III  and  IV  are  fairly  comparable  in  terms  of  solution  quality.  Recall, 
however,  the  eight  problems  for  which  Code  IV  utilized  excessive  CPU  time.  A 
user  of  Code  IV  is  advised  to  impose  a limit  on  execution  times  to  avoid  this 
occurrence.  Furthermore,  to  insure  replicability  when  several  problems  are 
solved  sequentially  by  Code  IV,  the  user  is  advised  to  reinitialize  the  random 
number  generator  between  problems  as  discussed  in  Section  II.  We  note  in 
passing  that  the  use  of  a double-precision  variable  in  Code  III  to  force  the 
double-precision  accumulation  of  inner  products  appears  to  be  a cost-effective 
approach  to  improving  solution  quality. 

Based  on  the  observation  that  interpolation  points  are  nearly  identical  for 
different  functional  forms,  we  anticipate  that  additional  investigation  of 
this  phenomenon  may  well  lead  to  a recommended  "starting  basis"  for  this  type 
of  problem.  The  theoretical  explanation  for  the  continuous  case  suggests  the 
use  of  the  roots  of  Chebyshev  polynomials  as  a starting  solution.  Further 
study  should  be  conducted  to  determine  the  practical  implications  of  this  sug- 
gestion in  the  discrete  case. 

We  caution  the  reader  that  the  results  obtained  from  this  experiment  are  valid 
only  for  the  type  of  problem  described  herein.  Further  evaluation  of  code 
performance  on  polynomial  approximation  problems  should  address  the  situation 
in  which  it  is  desired  to  find  the  best  polynomial  approximation  to  observed 
data.  This  type  of  problem  can  be  generated  by  P0LY2  and  should  be  utilized 
in  future  code  comparison  efforts.  The  importance  of  examining  many  class- 
es of  problems  is  emphasized  by  the  difference  in  the  results  reported  here 
and  those  in  earlier  experiments  [14]. 
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The  evaluation  described  in  this  paper  was  conducted  in  accordance  with  a sys- 
tematic methodology  for  obtaining  statistically  valid  results.  The  set  of 
test  problems  was  randomly  generated  within  a well-defined,  controlled  struc- 
ture known  to  arise  in  real  world  applications.  Performance  measures  were 
identified  for  both  solution  quality  and  computational  efficiency.  Finally, 
the  conditions  under  which  CPU  time  was  measured  were  chosen  to  minimize  vari- 
ability. The  results  of  the  experiment,  therefore,  permit  a valid  ranking 
among  the  four  codes  with  respect  to  code  performance  on  the  class  of  problems 
covered  by  this  experiment.  Although  the  rankings  may  vary  in  another  comput- 
ing environment,  the  experimental  approach  can  be  replicated  to  obtain  equally 
valid  results.  The  general  type  of  methodology  presented  in  this  paper  should 
serve  as  a guide  for  future  code  comparison  experiments. 
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